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In  this  report  we  consider  the  solution  to  integral  equations  that  arise  in  the  theory 
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operator  is  determined.  This  may  be  used  to  obtain  the  estimated  solution  to  an 
arbitrary  number  of  excitations.  It  is  shown  chat  the  outer  product  of  the  vectors 
needs  to  be  computed  to  compute  the  inverse  operator.  Since  this  is  an  expensive 
operation,  an  alternative  approach  is  determined  wherey  solution  to  multiple 
excitations  is  determined  without  computing  the  outer  products.  This  theory,  developed 
in  N-dimensional  space,  is  then  generalized  to  infinite  dimensional  space.  This  leads 
Co  a  definition  and  a  computational  procedure  for  inverse  integral  operator.  The 
theory  is  illustrated  by  computing  the  induced  currents  on  a  square  cylinder  for 
several  excitations. 
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SECTION  1 
INTRODUCTION 

The  techniques  that  are  available  to  compute  scattered  fields  from 
arbitrary  targets  may  be  classified  as  a)  analytical  techniques  and  b)  numerical 
techniques.  The  analytical  techniques  are  largely  based  on  asymptotic  theory 
ami  are  valid  when  the  target  size  is  sufficiently  large  compared  to  a 
wavelength.  These  techniques  are  limited  to  relatively  simple  targets  and  are 
limited  to  situations  where  the  requisite  diffraction  coefficients  are  available. 

The  numerical  techniques,  based  on  extensive  use  of  the  digital  computer, 
have  a  wide  applicability  and  the  target  configuration  is  arbitrary.  These 
techniques  formulate  the  problem  either  as  an  integral  equation  or  differential 
equation  subject  to  boundary  conditions  to  be  solved.  We  consider  in  this 
report,  the  solution  of  integral  equations. 

The  integral  equations  may  be  solved  directly  as  is  done  in  [1]  without 
explicitly  reducing  it  to  a  matrix  equation.  However,  with  this  method,  each 
time  the  excitation  is  changed,  the  solution  process  must  be  started  all  over 
again.  Thus,  when  there  are  a  large  number  of  excitations  to  be  considered,  as 
is  often  the  case,  presently  available  direct  solutions  are  inefficient.  A  more 
common  approach  is  to  take  a  projection  of  the  integral  equation  into  a  finite 
dimensional  space  reducing  it  to  a  matrix  equation  [2].  This  matrix  equation 
may  be  solved  by  several  conventional  methods  (3).  An  advantage  of  this  method 
is  the  availability  of  the  inverse  operator.  Once  computed,  this  may  be  used 
repeatedly  on  any  ntunber  of  excitations  simply  and  effectively.  This  approach 
is  clearly  among  the  best  to  solve  problems  with  multiple  excitations. 
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As  the  electrical  size  of  the  target  becomes  large,  the  order  of  the  matrix 
becomes  large  and  obtaining  the  inverse  operator  becomes  numerically  in¬ 
efficient.  In  addition,  round-off  error  tends  to  vitiate  the  computations  as 
the  order  of  the  matrix  becomes  large.  Iterative  techniques  have  been  used  to 
overcome  this  problem  [7]  -  [11].  However,  iterative  methods  whether  based  on 
Lanczos  iteration  [12]  or  conjugate  gradients  [13]  are  limited  to  only  one 
excitation.  If  the  excitation  changes,  the  problem  must  be  solved  in  its 
entirety  again. 

In  this  report  we  set  ourselves  the  task  of  devising  an  iterative  scheme 
that  may  be  used  with  not  just  one  excitation  but  an  arbitrary  number  of 
excitations.  That  is,  in  effect,  we  seek  to  compute  an  inverse  operator 
iteratively . 

Let  us  consider  a  solution  of  the  matrix  equation. 

Ax  -  F  (1.1) 

where  A  is  an  NxN  matrix.  Let  A  be  self  adjoint.  The  matrix  A  may  be  considered 
to  be  a  collection  of  N-Vectors,  j‘’*'  vector  being  simply  the  j*^^  column  of  the 
matrix.  Equation  (1.1)  may  be  put  in  the  form, 

<  Aj,  X  >  -  Vj  ,  j  =  1,2  N  (1.2) 

A  geometrical  interpretation  of  equation  (1.2)  is  to  find  the  vector  x  such  that 
its  projections  along  the  columns  of  the  matrix  (  <  Aj ,  x  >  )  equal  the 
components  of  the  excitation.  Th,is  problem  is  best  solved  by  setting  up  an 
orthogonal  co-ordinate  system.  This  is  accomplished  by  starting  from  an 
arbitrary  initial  vector  x^,  and  computing  a  sequence  of  linearly  independent  set 
of  vectors,  x^,  Ax^ ,  A^x^  ....  a”  x^.  Ihis  set  of  vectors  is  called  the  Krylov 
sequence  of  vectors  and  each  member  of  this  sequence  is  computed  by  operating 
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on  the  preceding  vector  by  the  matrix  A.  That  is,  Krylov  sequence  may  be 
computed  iteratively. 

The  Krylov  sequence  is  not  orthogonal.  These  vectors  may  be  orthogonalized 
by  the  Gram-schmidt  procedure  [3].  However,  this  procedure  is  numerically 
inefficient.  Lanczos  [12],  [6]  had  shown  that  it  is  only  necessary  to  enforce 
orthogonality  of  a  vector  from  Krylov  sequence  to  the  two  preceding  orthogonal 
vectors;  orthogonality  to  the  rest  of  the  vectors  is  assured.  Thus,  a  set  of 
A-orthogonal  vectors  are  computed  iteratively,  satisfying  the  relation. 

<  qi,  Aqj  >  -  <  Aqj,  qj  >  -  Ai  .  (1.3) 

It  is  shown  in  this  report  that  the  inverse  operator  may  be  thought  of  as 
consisting  of  N-components .  Each  component  matrix  is  shown  to  be  related  to  one 
of  the  A-orthogonal  vectors.  For  instance,  the  j*-*’  orthongonal  vector  qj  leads 
to  a  matrix  given  by,  Aj'^  q^  q^**.  Thus,  as  each  new  A-orthogonal  vector  becomes 
available,  the  inverse  operator  may  be  up-dated.  This  leads  to  an  interative 
computation  of  the  inverse  operator.  However,  these  computations  involve 
computing  the  outerproduct  qj  qj**.  Such  a  computation  is  numerically  expensive. 
A  technique  is  shown  whereby  solution  to  multiple  excitations  is  obtained  but 
without  explicitly  computing  the  outer  product. 

The  theory  developed  for  the  finite  dimensional  operator  is  then  extended 
to  infinite  dimensional  space.  This  leads  to  the  definition  of  the  inverse 
integral  operator.  This  inverse  operator  is  related  to  a  set  of  orthogonal 
functions  in  the  same  way  the  matrix  inverse  operator  is  related  to  the  q- 
vectors . 

In  section- 2  we  present  the  detailed  development  of  the  theory  of  the 
inverse  operator.  Section- 3  presents  a  numerical  implementation  of  the  theory 
and  a  few  examples.  Section-A  contains  a  discussion  of  the  results  and  some 
suggestions  for  further  research. 
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SECTION  2 


THEORY 

Here  we  consider  the  solution  of  Integral  Equations  (IE) .  In  the  field 
of  electromagnetic  scattering  Fredholm's  integral  equations  of  both  the  first 
and  second  kind  arise.  The  electric  field  integral  equation  (EFIE)  is  an 
integral  equation  of  the  first  kind  while  the  magnetic  field  integral  equation 
(MFIE)  is  of  the  second  kind.  The  general  form  of  these  IE's  is  given  below: 

g  -  Kf  (2.1) 

f  -  g  +  Kf  .  (2.2) 

I 

K  is  a  linear  integral  operator.  Equation  (2.2)  may  be  rearranged  so  that  it 
assumes  the  form  of  an  integral  equation  of  the  first  kind: 

g  -  f  -  Kf 
=  (I-K)  f 

-  K'  f  .  (2.3) 

Where  K’  is  a  new  operator  given  by 

K'  =  I  -  K 

Equation  (2.3)  is  an  integral  equation  of  the  first  kind.  The  kernel  of  this 
new  IE  involves  a  dirac  delta  function.  Thus,  without  loss  of  generality,  we 
consider  only  operator  equations  of  the  first  kind. 

The  integral  equation  may  be  thought  of  as  a  problem  in  infinite 
dimensional  space.  By  using  an  appropriate  technique  such  as  the  method  of 

moments  [2],  it  is  possible  to  take  a  projection  of  the  IE  (eq.  2-1)  frcir 

infinite  dimensional  to  an  N-dimensional  space.  This  procedure  reduces  the  " 
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to  a  finite  dimensional  matrix  operator  equation,  of  the  form  shown  below. 

Ax  -  F  (2.4) 

Where  the  operator  A  is  an  NxN  matrix. 

We  will  develop  the  theory  in  terms  of  the  finite  dimensional  operator  A 
and  eventually  obtain  the  solution  to  the  IE  as  the  limit  of  the  solution  to 
equation  (2.4)  as  N  -*  ®  .  The  infinite  dimensional  space  is  considered  to  be 
a  limiting  case  of  N-dimensional  space. 

Let  X  and  Y  be  two  vectors  defined  by 


X  -  [  Xi,  X2  . x„l' 

and  Y  -  (  y^,  y2  . y^]  ' 


The  prime  is  a  transpose  operator.  Then,  the  inner  product  of  two  vectors  is 
defined  by, 

<  X,Y  >  -  X"*  Y 

The  superscript  H  indicates  Hermitian  operation  and  X”  is  obtained  from  X  bv 
conjugating  and  transposing  the  elements  of  X. 

N  . 

<  X,Y  >  -  X«  Y  -  I  Xk  Yk  (2.5) 

k-1 

* 

where  Xj.  is  the  complex  conjugate  of  x^.  The  outer  product  of  two  vectors,  x  and 
y  is  an  NxN  matrix  Z,  given  by 

Z  -  X  Y“  .  (2.6) 

An  operator  is  said  to  be  self-adjoint  if  the  following  identity  holds. 

<  AX,  Y  >  -  <  X,  AY  >  (2.7) 
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Without  loss  of  generality  we  assume  the  operators  K  and  A  to  be  self-adjoint. 
The  impedance  matrix  that  is  obtained  in  electromagnetic  problems  is  not  self 
adjoint.  However,  the  problem  may  be  reformulated  to  yield  an  operator  that  is 
self  adjoint. 

ZI  -  V  (2.8) 

Multiplying  both  sides  of  equation  (2.7)  by  Z^,  we  obtain, 

(  Z«  Z  )  I  =  Z”  V  .  (2.9) 

It  may  easily  be  seen  that  the  operator  Z^  Z  is  self  adjoint. 

2.1  A- ORTHOGONAL  VECTORS 

We  now  determine  a  sequence  of  vectors,  q^,  q2,  q2  •  .  •  .  Qn  starting  from  an 
arbitrary  initial  vector  x^.  These  vectors  are  A-orthogonal  vectors.  That  is, 

<  Aqi,  qj  >  -  Ai  .  (2.10) 

Since  the  operator  A  is  self-adjoint, 

<  Aq^,  qj,  >  -  <  q^,  Aqj,  >  -  Ai  .  (2.11) 

These  vectors  are  obtained  by  A-orthogonalizing  the  Krylov  sequence  of  vectors. 
Starting  from  an  initial  arbitrary  vector,  x^,  the  Krylov  sequence  is  obtained 

as  Xq,  AXq  A^Xg  .  ,  A”'^  Xg.  Note  that  each  vector  of  this  sequence  is 

obtained  from  the  preceding  vector  by  operating  on  it  once  with  the  operator  A. 

The  set  of  Krylov  vectors  are  not  A-orthogonal.  These  may  be 
orthogonalized  by  a  special  procedure  such  as  the  Gram-Schmidt  procedure. 
However,  Gram-Schmidt  procedure  is  numerically  very  expensive.  Furthermore, 
due  to  the  special  properties  of  the  Krylov  sequence  of  vectors,  it  may  be  shown 
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that  orthogonal  vector,  qp,  may  be  obtained  by  ensuring  that  is 
orthogonal  to  only  the  two  preceding  vectors,  qp.j^  and  qp.2!  (Gram- Schmidt  would 
have  required  that  be  orthogonal  to  all  the  preceding  vectors.)  This 
clearly  results  in  an  efficient  computational  procedure.  We  now  give  the  details 
of  construction  of  the  q-vectors. 


qi  =•  Xo 

^2  “  ■  7  qi 


<  Aqi,  Aqj  > 

7i  “  - 

<  q^,  Aqj  > 


(2.12) 


qin  “  Aq^  -  7i  qi  -  Si  q,-i  .  i  -  2.3  .  N-1  (2.13) 


<  Aq^, 

Aqi  > 

<  qi- 

Aqi  > 

<  Aq^, 

Aqi-i  > 

<  qi-i. 

qi-i  > 

The  vectors  so  computed  constitute  a  complete  A-orthogonal  set  as  proved  in 
[1].  Note  that  the  computational  process  may  be  so  arranged  that  each  of 
these  vectors  are  computed  iteratively,  one  per  iteration.  The  A- 
orthogonality  of  these  vectors  may  be  expressed,  compactly,  using  matrix 
notation  as  follows: 

A  Q  -  A  (2.14) 


Where  is  the  transpose  of  the  NxN  matrix  Q,  Q  is  a  matrix  with  the  q-vectors 
as  its  elements.  That  is, 


Q  -  [  qi-  qa.  qs  .  qn  1 


(2.15) 
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A  is  a  diagonal  matrix,  whose  diagonal  elements  are  given  by 


At  -  <  qi  .  Aqi  > 


(2.16) 


2 . 2  INVERSE  OPERATOR 

We  now  obtain  explicit  expressions  for  the  inverse  of  the  operator  ir 
terms  of  the  q-vectors.  The  matrices  Q  and  Q®  are  partitioned  as  follows. 

Q  -  [  Qp  Qp  ]  .  (2.17) 

Where  Qp  is  an  N  x  P  matrix  obtained  from  the  first  P  A- orthogonal  vectors. 
That  is, 


Qp  -  [  qi  q2  • • •  qp  ]  •  (2.i8) 

Qp  is  an  N  X  (N-P)  matrix  defined  below. 

Qp  -  (  Vi  V2  •  •  •  %  ]  (2.19) 

Qp'*  and  Q^®  are  matrices  obtained  from  transposing  and  conjugating  Qp  and  Qp. 
From  these  matrices,  Q®  may  be  shown  to  be, 


LQp“  J 

Now,  taking  the  inverse  of  equation  (2.1A)  gives 

[  qH  A  Q  ]-'  -  A-' 

Q'^  A'^  [Q“]'^  -  A'^ 


(2.20) 


(2.21) 


A'^  -  Q  A'^  Q® 


(2.22) 
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Note  that  since  A  is  a  diagonal  matrix,  A*^  is  simply  another  diagonal  matrix, 
whose  diagonal  elements  are  Aj^*^.  Thus,  equation  (2.22)  shows  the  relationship 
between  the  inverse  operator  and  the  A-orthogonal  vectors.  However,  equation 
(2.22)  implies  that  the  inverse  operator  may  be  computed  only  after  all  the  q- 
vectors  are  computed.  In  fact,  the  q-vectors  are  computed  iteratively;  an 
additional  q-vector  is  available  after  each  iteration.  It  can  be  shown  that 
the  inverse  operator  itself  may  be  obtained  iteratively.  The  diagonal  matrix 
may  be  partitioned  as  shown  below. 


A-'  - 


(2.23; 


The  Ap’^  and  A^'^  are  diagonal  matrices  of  orders  PxP  and  P  x  P  (P  -  N  -  P)  . 
Then,  from  equation  (2.22), 


A-'  -  [  Qp  Qj  ] 


-i 


0 


■  ■ 

.  Qp“  . 

(2.24) 


-  (  Qp  Qp  ] 


-1  n  H 
p 


Ap-^  Q, 


A-'^  0-® 

“p  ^P 


[  Qp  Ap-'  Qp«  +  Qp«  QpS  ] 

[  Bp  +  Bp.  1 


(2.25) 


Where , 


and 


Bo  -  [  Qp  Ap-‘  Qp”  ] 


Bp  -  [  Qp  Qp«  ] 


(2.26) 

(2.27) 


Bp  may  be  thought  of  as  the  best  estimate  of  the  inverse  operator  after  P 
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iterations . 


Indeed,  each  A-orthogonal  vector  q  corresponds  to  a  matrix 


component  of  the  inverse  operator.  By  direct  expansion  of  equation  (2.26), 

Bp  -  [qi  qi”  +  q2  A2'^  qa"  +  -  +  Qp  Ap'^  qp'’] 

P 

-  2  Aj-'  qj  q/  .  (2.28) 

Thus,  the  outerproduct  (a  matrix)  formed  from  the  vector  q^  is  a  component 

of  the  inverse  operator.  It  must  be  pointed  out  that  computing  an 
outerproduct  is  numerically  a  very  expensive  operation.  Indeed,  solution  to 
the  original  operator  problem  for  multiple  inputs  may  be  computed  without 
actually  computing  the  outerproduct  and  hence  the  inverse  operator! 


2.3  SOLUTION  WITHOUT  EXPLICIT  INVERSE  OPERATOR 


We  now  consider  solution  to  the  operator  problem  with  M-excitations . 


AXj  -  pj  j  -  1,2,  _ M  .  (2.29) 


Let  P  iterations  of  the  procedure  be  carried  out  resulting  in  the  computation 

j 

of  P  A-orthogonal  vectors  q^,  q2,  ....  qp.  Then  the  solution  Xp  corresponding 
to  the  excitation  F-^  is  given  by 


P  l 

-  [  s  —  qt  qi”  ]  pj 

i-1  A* 


P  1 

2  —  qt  qt”  Fj 

i-l  A* 


1 

-  <  qt ,  pj  >  q(  (2.30) 

A* 


P 

S 

i-1 
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p  <  > 

-  2  -  qt  (2.31) 

■*“1  <  qt.  Aq^  > 

From  equation  (2.31)  it  is  evident  that  it  is  only  necessary  to  compute  the 
inner  products  <  q^,  Fj  >  and  not  the  outerproduct  q^  q^®  to  obtain  the 
solution.  This  is  clearly  much  more  efficient  than  computing  the  inverse 
operator  directly.  Even  for  purposes  of  storing  the  inverse,  it  is  best  to 
store  the  collection  of  q-vectors,  the  Q-matrix,  rather  than  A  itself. 
Equation  (2.31)  may  also  be  written  as, 

<  q«.  Fi  >  <  qp,  Fj  > 

-  qt  +  -  qp 

<  q«.  Aq^  >  <  qp,  Aqp  > 

<  qp.  > 

-  Vi  +  - 

<  qp.  Aqp  >  (2.32) 

Thus,  after  P  iterations,  as  the  vector  becomes  available,  the  solution 
computed  at  the  end  of  P-1  iterations  is  up-dated  by  the  addition  of  the 
factor  <  qp,  >  qp  /  <  qp,  Aqp  >.  This  is  analogous  to  up-dating  the  inverse 
operator  by  the  addition  of  the  matrix  Ap'^^  qp  qp®. 

2 . 4  ITERATIVE  SCHEME 

With  the  theory  in  hand,  we  are  now  in  a  position  to  outline  the 
iterative  scheme  to  be  implemented  on  a  digital  computer.  Let  x,,  be  an 
arbitrarily  chosen  initial  vector. 

Step  1; 

qi  -  Xo 

Ai  -  <  q^,  Aq^  > 


p-1 
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1 

-  -  <  qi.  F-i  >  ,  j  -  1,2,  ...  M 

,  J  -  1.2,  ...  H  (2.33) 

Step  2: 

1 

7i  -  -  <  Aqi  ,  Aq^  > 

Ai 

qa  ”  Aqi  -  7i  qi 

A2  “  <  qa ,  Aq2  > 

1 

ta'’  - - <  qa .  FJ  >  q2  .  j  -  1 . 2  .  .  .  M 

Aa 

X2^  -  Xi-i  +  t2j  ,  j  -  1.2  ...  M  .  (2.34) 


Step  3:  For  i  -  2,  3  . 

1 

7i  =  -  <  Aqj,  Aq^  > 

Ai 


1 


- 

-  <  Aq^, 

Ai-i 

Aqi 

qm  ” 

Aqi  -  7i  qi  - 

Aj+i  ” 

<  qiti-  Aq^^i 

> 

j 

tin  -  -  <  qin.  FJ  >  qi^i  ,  j  -  1.2.  ...  M 


Xin  -  Xi-j  +  ti^i'^  ,  j  -  1,2  ...  M 


e-j  -  llpj  -  Ax^^ill  .  j  -  1,2  ...  M  .  (2.35) 
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Step  4: 


Terminate  if  is  sufficiently  small  or  go  back  to  Step- 3. 

This  procedure  computes  one  A-orthogonal  vector  for  each  iteration  and  uses 
this  vector  to  update  this  solution  for  all  M  excitations  and  not  just  one  excitation. 
Note  that  the  inverse  operator  A'^  is  not  computed  explicitly.  As  observed  earlier, 
the  matrix  Q,  which  is  a  collection  of  the  A-orthogonal  vectors  serves  the  same 
purpose  as  A'^  but  numerically  is  much  more  efficient. 

2.5  INVERSE  INTEGRAL  OPERATOR 

The  theory  of  the  inverse  operator  has  been  developed  in  N-dimensional  space  . 
We  now  let  the  number  of  dimensions,  N,  approach  infinity  so  that  the  matrix  operator 
becomes  integral  operator.  It  is  worth  remembering  that  the  matrix  operator  which 
becomes  integral  operator  was  originally  obtained  as  a  projection  of  the  integral 
operator  into  finite  dimensional  space.  That  is, 

asN-»<«>,A-*K  (2.36) 

The  vectors  of  finite  dimensional  space  become  functions  in  infinite  dimensional 
space.  That  is, 

qi  -  (2.37) 

The  inner  product  in  infinite  dimensional  space  of  two  functions  f  and  g  may  be 
defined  as, 

<  f,  g  >  -  /  *f  g  di  (2.38) 

As  before,  the  operator  K  is  self  adjoint  if, 

<  K  f,  g  >  -  <  f,  Kg  >  (2.39) 

We  now  consider  the  solution  of  the  integral  operator  equation,  for  multiple 
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excitations : 


K  f j  -  g->  .  j  -  1,2  ...  M  .  (2.40) 

As  in  the  case  of  finite  dimensional  space,  the  functions  are  computed  such 
that  they  satisfy  the  orthogonality  condition. 


<  Ci  .  >  -  Ai  Sjj  (2.41) 

The  difference  is  that,  in  N-dimensional  space  the  number  of  q- vectors  are  finite 
in  number,  equal  to  N  but  in  infinite  dimensional  space,  are  infinite  in  number. 
However,  the  inverse  operator  K'^  is  related  to  the  in  a  manner  analogous  to 
the  way  A‘^  is  related  to  q^. 

In  N-dimensional  space,  the  solution  after  P- iterations  is  given  by  equation 

(2.31): 

p  <  qj,  F-*  > 

XpJ  -  Bp  Pj  -  S  -  qe 

i-1  <  qe.  Aqe  > 

Generalizing  this  equation  to  infinite  dimensional  space. 


P  <  > 

V  ■  s'*  -  2  -  Ct  .  j  -  1-2,  .  .  .  M  (2.42) 

f-1  <  Ce  > 


Thus,  obtaining  an  iterative  solution  to  integral  equation  with  multiple  excitations 
consists  of  the  following  steps. 

Step  1: 

Determine  K-orthogonal  functins  starting  from  an  arbitrary  initial 
function  fo  and  Krylov  sequence  of  functions. 

Step  2: 

These  are  computed  iteratively,  ie,  one  function  during  each 
iteration.  This  function  may  be  used  to  compute  the  update  for  the 
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solution  as  given  by  equation  (2.42).  The  inverse  operator  may  be 
expressed  as  follows: 


“  it  it^ 

K"^  -  S  - 

f-1  <  KCt  >  (2.43) 

The  outerproduct  of  two  functions  results  in  a  new  operator.  The  operation 

of  this  operator  on  a  function  f  is  defined  by: 

(  ft  )  f  -  <  ft.  f  > 

-  <  ft.  f  >  ft  (2.44) 

The  effect  of  this  operation  is  to  transform  the  function  f  into  a  scalar  multiple 
of 

As  may  be  seen  from  equation  (2.43),  this  inverse  operator  depends  only  on 
^'s  and  the  operator  K  and  is  independent  of  excitation.  We  now  give  an  iterative 
scheme  for  the  solution  of  equation  (2.40). 

2.5.1  ITERATIVE  SCHEME  FOR  INVERSE  INTEGRAL  OPERATOR 

The  computational  scheme  for  the  solution  of  the  integral  equation  with 
multiple  right  hand  sides  is  given  here. 

Step  1.  Choose  an  arbitrary  initial  function  f^.  Then, 

$i-fo  ; 

1 

- <  fi.  g-*  >  f  .  j  -  1,2  ...  M  (2.45) 

Ai 

Step  2. 

1 

7i  -  —  <  K^i.  > 

Ai 
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^2  -  -  7i 


Az  -  <  ^2.  K?2  > 
1 


£2^  -  +  -  <  ?2  .  g-^  >  ^2 

A2 

Step  3. 


7i  - 

1 

—  <  KCi  . 

At 

> 

= 

1 

—  <  . 
Ai-1 

K^i-i  > 

“  Ui  -  7i 

- 

Ai+i 

J 

^i+l 

1 

-  fi^  +  - 

Ai+i 

gj  > 

ej  - 

II  -  Kfinll 

t 

j  “ 

j  -  1,2,  . . .  M  .  (2.46) 


?  i+i  .  j"*l,2,...M 

1,2.  . . .  M  .  (2.47) 


Step  4. 

If  e-*  is  not  sufficiently  small,  go  back  to  Step-3. 

It  may  be  noted  that  in  order  to  obtain  the  solution  to  the  multiple  excitation 
problem,  it  is  only  necessary  to  store  at  most  three  ^'s  at  any  given  time.  However, 
if  the  inverse  operator  were  to  be  stored  to  be  used  at  a  later  time  on  more 
excitations,  then  all  the  ?'s  computed  need  to  be  stored. 

2 . 6  DEGENERACY 

The  iterative  scheme  described  in  Sections  2.3  and  2.5.1  may  terminate 
prematurely  without  the  error  becoming  sufficiently  low  if  the  initial  vector  x. 
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or  fg  are  deficient.  The  initial  vector  must  contain  components  along  all  the 
eigenvectors  of  the  operator  A.  Otherwise  the  Krylov  sequence  ceases  to  produce 
new  linearly  independent  vectors  after  a  certain  number  of  iterations.  When  this 
happens ,  it  is  not  possible  to  compute  a  new  A-orthogonal  vector  leading  to  a  premature 
termination  of  the  iterative  scheme.  Similar  comments  apply  for  the  case  of  continuous 
operator  as  well.  In  such  cases  of  degeneracy,  the  iteration  may  be  re-started 
using  the  following  procedure. 

Let  be  the  initial  vector.  Let  P  A-orthogonal  vectors  q^,  q2  ...  Qp  be 
computed  before  degeneracy  occurred.  Then  the  iterative  procedure  may  be  continued 
from  another  vector  x^^.  However,  x^^  must  be  A-orthogonal  to  all  the  preceding 
P  vectors.  This  may  be  ensured  in  the  following  way.  Define  qp+j  such  that, 

P  1 

Vi  -  -  2  —  <  x„2,  Aqt  >  qZ  (2.48) 

i-1  At 

“Ip+i  nisy  be  shown  to  be  A-orthogonal  to  all  the  preceding  q-vectors.  That  is, 

<  Vi-  Aqj>-0,  j-l,2...P  .  (2.49) 


P  1 

<  qp+i.  Aqj  >  -  <  x„2,  Aqj  >  -  Z  —  <  ,  Aq^  >  <  q^,  Aq,  > 

i-1  Aj 


-  <  x„2,Aqj  >  -  Z  —  <  x„2,  Aq^  >  Aj 
i-1  At 


With  qp^.^  defined  this  way,  the  iteration  may  be  continued  from  Step-3  of  Section 
2.4.  However,  it  must  be  noted  that  equation  (2.48)  implies  that  P  A-orthogonal 
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vectors  are  stored  and  are  available.  It  is  possible  to  devise  a  scheme  wherebv 
it  is  not  necessary  to  store  the  A-orthogonal  vectors.  Two  vectors  and 
are  initially  chosen  and  the  iteration  is  started  with  x^^  as  given  in  section 
2.4.  In  addition,  however,  after  each  iteration,  x^^  is  modified  so  that  it  is 
A-orthogonal  to  q-vector  that  is  computed.  During  iteration,  then,  x^^  is 
replaced  by, 

1 

Xq^  -*  x„2  -  —  <  qt  >  q^  (2.51) 

At 

This  procedure  ensures  that  x^^  is  always  A-orthogonal  to  all  the  available  q's. 
Of  course,  similar  procedure  may  be  used  for  the  continuous  operator  case  as  well. 


SECTION  3 


NUMERICAL  ANALYSIS 

The  numerical  Implementation  of  the  theory  described  in  Section- 2  is  discussed 
here.  We  present  here  the  numerical  considerations  that  result  in  the  most  efficient 
code  as  well  as  the  numerical  difficulties  that  may  arise  and  some  prescriptions 
to  overcome  these  problems.  Some  numerical  examples  are  presented. 

3.1  NUMERICAL  CONSIDERATIONS 

a)  A  significant  part  of  the  computation  consists  in  computing  from 

Aq^  by  subtracting  the  projections  of  q^  and  qj-i.  This  may  be  accomplished  by 
computing  and  5^  and  computing, 

qin  -  Aq^  -  7i  qi  -  (3.1) 

as  described  in  Section  2.6.1,  Step-3.  A  better  way  to  accomplish  this  is  to  carry 
out  the  recursive  procedure  described  below  [3],  [4]. 

qi-i  -  Aqi  -  7i  qi  (3.2) 

7i  “  <  Aq^ ,  Aqi  >  /  Ai  (3.3) 

qi+i  ”  qi+i  ■  qi-i  (3.4) 

-  <  qi+i.  Aqi-i  >  /  Ai  (3.5) 

The  projections  7^  q^  and  6^  q^.^  are  removed  one  at  a  time.  In  a  finite  precision 
arthmetric,  this  leads  to  a  more  stable  computation. 

b)  The  impedance  matrix  Z  that  comes  up  in  electromagnetic  problems  is  not 
self  adjoint.  However,  the  problem  may  be  reformulated  to  yield  a  new  matrix 
operator  that  is  self  adjoint.  That  is,  'the  original  matrix  equation, 
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Z  I  -  V 


(3.6) 


is  pre-multiplied  by  the  matrix  7^  to  obtain, 

(  Z»  Z  )  I  -  z“  V  .  (3.7) 

Equation  (3.7)  may  be  put  in  the  form 


Ax  -  F 

Where , 

A  -  Z*^  Z 


(3.8) 


(3.9) 


Multiplication  of  a  matrix  by  another  matrix  is  numerically  a  very  expensive 
operation.  In  fact,  the  numerical  scheme  may  be  arranged  so  that  one  never 
needs  to  carry  out  the  operation  in  equation  (3.9).  The  entire  iterative 
scheme  discussed  in  section  2  is  based  on  a  matrix  operating  on  a  vector, 
AX.  The  operation  (Z*'  Z)  x  may  be  carried  out  in  two  steps  as  follows. 
Let , 


Y  -  (  Z«  Z  )  X 

(3, 

,10) 

then. 

-  Zx 

(3, 

.11) 

and 

Y  -  Z“  Yi 

(3. 

.12) 

That  is,  the  result  required  in  equation  (3.10)  is  obtained  by  implementing 
equations  (3.11)  and  (3.12).  Both  of  these  equations  consist  of  matrix 
operating  on  a  vector,  a  numerically  efficient  operation. 


c)  In  an  infinite  precision  arithmetic,  the  q-vectors  computed  using  the 
procedure  in  section  2.5  will  always  be  A-orthogonal .  However,  in  a  finite 
precision  arithmetic,  as  invariably  is  the  case,  the  q-vectors  slowly  lose 
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their  orthogonality.  Such  a  loss  of  orthogonality  completely  destroys  the 
validity  of  computations  and  thus  must  be  monitored  and  must  be  prevented 
from  happening.  In  this  investigation,  we  monitor  the  loss  of  global  or¬ 
thogonality  by  computing  a  constant  "C"  as  defined  below,  after  computing 

*?p+i  • 

C  -  <  Aqi  ,  qp^i  >  (3.14) 

If  C  were  to  be  more  than  a  prescribed  value,  e,  (say  0.05)  ,  qp+i  is  deemed 
to  have  lost  orthogonality  and  is  rejected;  The  iteration  is  restarted  from  a  new 
vector  after  ensuring  that  it  is  orthogonal  to  all  the  preceding  P  number 
of  q-vectors. 

3.2  NUMERICAL  ILLUSTRATIONS 

The  target  is  a  square  cylinder  illuminated  by  a  TE  plane  wave.  The  angle 
of  incidence  is  arbitrary.  Each  angle  of  incidence  corresponds  to  a  different 
excitation.  For  purposes  of  illustration,  we  consider  three  angles  of  incidence 
ranging  from  90“  to  270“  as  shown  in  Figure-1.  We  consider  three  different  sizes 
with  W=0 . 2A  ,  0 . 3A  and  0 . 4A  .  Each  of  these  cylinders  is  illuminated  with  4>i  =  180“ , 
225“,  270“. 

A  standard  method  of  moments  code  [14]  is  used  to  generate  the  impedance 
matrix  Z  and  the  excitatin  vectors  VJ ,  j  -1,2,3.  The  matrix  equation  is  therefore  , 

Zl-v->  j-1,2,3  (3.15) 

This  equation  is  solved  using  the  procedure  described  in  this  report  iteratively 
and  for  all  the  excitations  simultaneously.  The  currents  obtained  using  this 
procedure  are  checked  by  comparing  them  to  currents  obtained  using  Grout's  method 

[5]. 

Figures  2-5  show  the  results  of  0.2A  square  cylinder.  Figures  2-4  show 
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that  the  currents  obtained  using  the  iterative  scheme  compare  well  with  the  currents 
obtained  using  elimination  procedure.  Figure  5  shows  the  monotonic  decrease  of 
the  r*m*s  error  as  iteration  progresses.  Similar  results  are  shown  in  figures 
6-9  for  0.3A  case  and  in  figures  10-13  for  0.4A  case. 

In  each  of  these  cases  the  iteration  is  started  by  choosing  the  initial  vector 

to  be 

x„  -  -  (3.16) 

<  AF^  F^  > 

F^  -  Z”  (3.17) 

The  vector  F^  is  the  right  hand  side  of  the  transformed  matrix  equation  corresponding 
to  —  180°.  Examination  of  figures  5,  9  and  13  show  that  error  decreases 
fastest  f  the  -  180°  case.  The  rate  of  convergence  seems  to  depend  on  the 
initial  choice  of  the  vector  and  the  best  choice  seems  to  be  the  right  hand  side 
of  the  matrix  equation  itself. 
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Figure  2.  Currents  on  a  perfectly  conducting  square  cylinder,  w=0.2X,  PHI=180 
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Figure  3.  Currents  on  a  perfectly  conducting  square  cylinder,  w*0.2X,  PHI=222. 


Figure  4.  Currents  on  a  perfectly  conducting  square  cylinder,  w=0.2A,  PHI=227 
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Figure  5.  RMS  Error  Vs  Iteration  Number  for  three  different  excitations  of  square  cylinder, 


Figure  7.  Currents  on  a  perfectly  conducting  square  cylinder,  w-O.JA,  PHI-222. 
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gure  9.  RMS  Error  Vs  Iteration  Number  for  three  different  excitations  of  square  cylinder,  W-0.3X. 
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Currents  on  a  perfectly  conducting  square  cylinder,  w=0.4X,  PHI=222. 
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SECTION  4 


DISCUSSION 

In  this  report  we  have  developed  the  theory  of  inverse  matrix  operator  and 
inverse  integral  operator.  We  have  shown  that  it  is  possible  to  compute  the 
solution  iteratively  for  multiple  excitations  rather  than  just  a  single 
excitation. 

While  the  theory  is  complete,  several  numerical  problems  exist.  The  basis 
of  the  inverse  operator  is  the  construction  of  a  proper  orthogonal  set  of 
vectors.  In  a  finite  precision  machine,  as  the  iteration  progresses, 
orthogonality  is  slowly  lost.  If  the  vectors  were  not  A-orthogonal ,  the  theory 
developed  is  not  applicable  and  the  compu’  tions  become  meaningless.  Hence, 
efforts  must  be  made  to  monitor  the  orthogonality  and  when  the  loss  of 
orthogonality  is  signaled,  the  iteration  must  be  re-started. 

Paige  [15]  had  shown  that  loss  of  orthogonality  is  not  simply  because  of 
accumulation  of  round-off  error.  Rather  it  is  because  of  the  combined  effects 
of  round-off  error  and  convergence  of  a  few  eigenvectors.  In  fact,  the  losses 
of  orthogonality  occurred  only  along  the  directions  of  converged  vectors.  We 
recommend  therefore  the  inclusion  of  selective  orthogonalization  rather  than 
complete  re-orthogonalization  which  is  numerically  expensive.  As  the  computation 
continues,  each  loss  of  orthogonality  becomes  a  signal  that  a  new  eigenvector 
is  available.  This  eigenvector  is  computed  and  is  removed  from  all  later  q's 
as  otherwise  an  image  of  the  eigenvector  appears  again  and  again. 

Much  further  work  needs  to  be  done  to  develop  the  inverse  operator  theory. 
It  would  be  very  interesting  indeed  to  compute  the  inverse  integral  operator; 
That  is,  implement  the  scheme  discussed  in  section  2.6.1.  If  this  were  to  be 
successful  it  completely  eliminates  the  need  to  explicitly  reduce  the  integral 
equation  to  matrix  equation. 
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ROME  LABORATORY 

Rome  Laboratory  plans  and  executes  an  interdisciplinary  program  in  re¬ 
search,  development,  test,  and  technology  transition  in  support  of  Air 

O 

Force  Command,  Control,  Communications  and  Intelligence  (C  l)  activities 
for  all  Air  Force  platforms.  It  also  executes  selected  acquisition  programs 
in  several  areas  of  expertise.  Technical  and  engineering  support  within 
areas  of  competence  is  provided  to  ESD  Program  Offices  (POs)  and  other 
ESD  elements  to  perform  effective  acquisition  of  C^I  systems.  In  addition, 
Rome  Laboratory's  technology  supports  other  AFSC  Product  Divisions,  the 
Air  Force  user  community,  and  other  DOD  and  non-DOD  agencies.  Rome 
Laboratory  maintains  technical  competence  and  research  programs  in  areas 
including,  but  not  limited  to,  communications,  command  and  control,  battle 
management,  intelligence  information  processing,  computational  sciences 
and  software  producibility,  wide  area  surveillance/sensors,  signal  proces¬ 
sing,  solid  state  sciences,  photonics,  electromagnetic  technology,  super¬ 
conductivity,  and  electronic  reliability/maintainability  and  testability. 


